Non-axisymmetric instability of shear-banded Taylor-Couette flow 
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Recent experiments show that shear-banded flows of semi-dilute worm-like micelles in Taylor- 
Couette geometry exhibit a flow instability in the form of Taylor-like vortices. Here we perform 
the non-axisymmetric linear stability analysis of the diffusive Johnson-Segalman model of shear 
banding and show that the nature of this instability depends on the applied shear rate. For the 
experimentally relevant parameters, we find that at the beginning of the stress plateau the instability 
is driven by the interface between the bands, while most of the stress plateau is occupied by the bulk 
instability of the high-shear-rate band. Our work significantly alters the recently proposed stability 
diagram of shear-banded flows based on axisymmetric analysis. 



PACS numbers: 83.80. Qr, 47.50.-d, 47.50. Gj 

Flows of complex fluids often differ dramatically from 
their Newtonian counterparts [1 j. One of the most strik- 
ing examples of such differences is the phenomenon of 
shear-banding which is widespread in flows of micellar so- 
lutions [U [2] , granular media [3 , foams [4] and colloidal 
glasses [5] . Its name is derived from the observation that 
when sheared these materials often split into two or more 
bands of different shear rates, viscosities and viscoelastic 
properties [TJ [2]. In semi-dilute solutions of worm-like 
micelles, shear-banding is often associated with the so- 
called shear-stress plateau - the range of applied shear 
rates j app for which the shear stress E is approximately 
independent of j app [2j [6] [7] . Outside the stress plateau, 
simple ID shear flow remains homogeneous, while for the 
values of j app on the plateau, it splits into two regions of 
simple shear flow with different shear rates [2]. 

Recent experiments with semi-dilute worm-like mi- 
cellar solutions have demonstrated that this picture of 
steady ID shear-banded flows only holds on average [8]- 
[TT] . Visualisation experiments in flows between two ro- 
tating coaxial cylinders, the Taylor-Couette flow, have 
shown that on the stress plateau, there exist large in- 
stantaneous fluctuations in local stresses and the local 
position of the interface between the two bands that are 
associated with a 3D hydrodynamic instability [8HTU]. 
This instability has the form of doughnut-shaped vor- 
tices threaded by the inner cylinder and stacked in the 
vorticity direction (Taylor-like vortices) and is mainly lo- 
calised in the high-shear-rate band [TT]. Its origin is not 
inertial due to the low flow velocities and large viscosities 
involved [2]. 

Presently, the precise nature of the instability driving 
the Taylor-like vortices and undulations of the interface 
remains unclear. Possible mechanisms include an insta- 
bility driven by the presence of the interface (interfa- 
cial mode) or a bulk instability inside one or both bands 
(bulk mode). The interfacial mode is a viscoelastic ana- 
logue of the Kelvin-Helmholtz instability and arises in 



co-flows of several liquids with a discontinuity in their 
normal-stress differences across the interface [12] . The 
bulk mode, which has been observed in Taylor-Couette 
flow of polymer solutions, is driven by large values of the 
first normal stress difference N± in the bulk, and the pres- 
ence of curved streamlines p~3j [14]. Unlike Newtonian 
Taylor-Couette flow where the first instability results in 
an axisymmetric stack of Taylor vortices, linear stabil- 
ity analysis predicts that viscoelastic Taylor-Couette flow 
is unstable towards non-axisymmetric (wavy in the az- 
imuthal direction) Taylor-like vortices [16]. 

Previous attempts to uncover the nature of the insta- 
bility in the shear-banded Taylor-Couette flow focused on 
the axisymmetric version of this flow [15] . By performing 
numerical simulations, Fielding has concluded that micel- 
lar solutions split into three groups depending on their 
material properties: fluids with low values of the first 
normal stress difference Ni in their high-shear-rate band 
exhibit the interfacial instability [16], while fluids with 
high values of Ni exhibit the bulk instability in curved 
geometries [15]. For intermediate values of iVi, the flow 
was predicted to be stable, and the precise position of 
this region was found to depend on the curvature of the 
Taylor-Couette cell [15]. The main problem with this sce- 
nario is that the fluid used by Lerouge et al. [8] [TO] [11] to 
observe the instability most probably belongs to the sta- 
ble region of the stability diagram proposed in [15 (see 
Supplemental Material), casting doubt on the proposed 
mechanism. 

In this Letter we show that this contradiction is re- 
solved if one relaxes the assumption of axisymmetric flow. 
We perform a linear stability analysis of the shear-banded 
Taylor-Couette flow with respect to non-axisymmetric 
disturbances for a model fluid with material properties 
similar to the solutions used in the experiments by Ler- 
ouge et al. [8] [10] [TT]. We find that while axisym- 
metric perturbations are stable as predicted by Fielding 
[T5] , non-axisymmetric modes are unstable for all applied 
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shear rates on the stress plateau. The nature of the in- 
stability depends on the applied shear rate: very close to 
the beginning of the plateau, we find this is dominated by 
the interfacial mode, while most of the plateau is unsta- 
ble via the bulk mode in agreement with the observations 
by Fardin et al. [TT] . 

Our model consists of the momentum conservation 
equation in the limit of negligible inertia, 



Vp + T] s Av + V • r = 0, 



(1) 



and the diffusive Johnson- Segalman (DJS) equation [TJ 
fT5HT8] for the viscoelastic stress r: 
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where r= dr/dt + t? • Vr — Vt^ • r — r • and 

r= dr/dt + v ■ Vr + • r + r • Vv^ are the upper- and 
lower-convected derivatives, respectively pQ; (Vv)ij = 
dvj/dxi is the velocity gradient tensor, and f denotes 
the transpose. Here, v and p are the velocity and pres- 
sure, T] s and rj p are the viscosities of the Newtonian "sol- 
vent" and viscoelastic (micellar) components, and A is 
the Maxwell relaxation time. The slip parameter a con- 
trols the degree of non-affine deformations under flow [T7] 
and Eq.Q predicts a non-monotonic flow curve (shear- 
banding) as long as a ^ ±1 and f] p /r]s > 8. The stress- 
diffusion term / 2 V 2 r in Eq.Q prevents stress variations 
on scales smaller than I and uniquely selects the value of 
the plateau shear stress [18]. The lengthscale I sets the 
size of structural correlations in the fluid, and controls 
the width of the interface between the bands and the 
extent of the near-wall layers where the bulk behaviour 
changes rapidly in order to match the boundary condi- 
tions. Finally, we assume that the fluid is incompressible, 
V • v = 0, and apply the no-slip boundary conditions for 
the velocity v and the no-flux boundary conditions for 
the stress tensor r [TO] [TO] . 

We consider the Taylor- Couette flow between two 
coaxial cylinders. The inner cylinder of radius R\ rotates 
with the angular velocity while the outer cylinder of 
radius R2 is kept stationary. The fluid flows in the gap of 
size d = R2 — R1 and the relative curvature of the Taylor- 
Couette cell is set by e = d/R\. We choose the parame- 
ters in our model, Eqs.( |l|2[ ), to closely match their exper- 
imental values as used by Lerouge et al. [8] [TO] [11] and 
set Ri = 13.33 mm and d = 1.13 mm. To find the relax- 
ation time A, the slip parameters a and the viscosities r] s 
and we fit the rheological predictions of the Johnson- 
Segalman model in plane homogeneous shear flow to the 
cone-and-plate measurements for a solution used in [TO] . 

In Figjl] we compare the experimental data (crosses, 
courtesy of S. Lerouge) against the fit to the Johnson- 
Segalman model (solid and dashed lines) with A = 0.51 8, 
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FIG. 1. Total shear stress £ and first normal-stress difference 
Ni as a function of the applied shear rate j apP - Crosses - 
experimental measurements by S. Lerouge et al. [10] in a 
cone-and-plate geometry; solid line - JS constitutive curve for 
planar homogeneous shear flow; dashed line - prediction of 
the JS model in planar shear-banded flow. 
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FIG. 2. (left) Maximum growth rate $te(a*) as a function of 
Wi on the stress plateau (boundaries denoted by the dashed 
lines) for various values of the diffusion length /: squares - / = 
73 jLtm, circles - / = 13 /xra, triangles -1 = 4 /am. (right) Axial 
and azimuthal wavelengths k* and m* of the most unstable 
mode as a function of Wi for / = 13 jim. The lines are drawn 
to guide the eye. 



a = 0.985, T] s = 1.1 Pa ■ s and r] p = 33.0 Pa ■ s. These 
parameters are similar to the linear rheological measure- 
ments of [10 , A « 0.23 s and r] s + r] p ~ 55 Pa • s. We 
mostly consider two values of the structural lengthscale, 
I = 13 fim and I = 4 /im, which are consistent with the 
experimentally determined values [TO] [TO] [20] . 

First we introduce cylindrical coordinates (r, z) and 
solve the equations of motion (T|2) numerically for the 
steady ID shear-banded profile using a pseudospectral 
Chebyshev-tau method [21] . To accurately resolve sharp 
interfaces we split the computational domain into three 
regions, two for each shear band and one corresponding 
to the interface, and use adaptive domain decomposition 
algorithm to discretize them independently. A similar 
method was employed in [22] . 

The linear stability analysis is performed by perturb- 
ing the ID profile with 3D perturbations of the form 
(v, r,p) ~ exp (at + ikz + zm#), linearising the equa- 
tions of motion, and solving the resulting generalised 
eigenvalue problem. Below, the applied shear rate is re- 
ported in terms of the dimensionless Weissenberg number 
Wi = Aft/(e(e + 2)). 
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FIG. 3. (top) Superimposed spectra showing the most unstable eigenvalue for the shear-banded (dots - I = 13 fim, triangles 
- I — 4 fim) and the auxiliary (crosses - I — 4 fim) setup, (a) Wi — 29, k = 20.54 and m — 21; (b) Wi — 3.0, k = 14.4 and 
m = 30. (bottom) Contour plots of the dispersion relation 9ftecr(fc, m) as a function of m and k. White triangle marks the most 
unstable eigenmode. (c) Wi — 29 and / = A fim, green dashed lines - contour plot for the auxiliary system; (d) Wi = 3.0 and 
/ = 13 jitm; (e) Wi = 6.4 and I = 4 fim showing two types of instability. 



In Figj2] we plot the real part of the most unstable 
eigenvalue a* and the corresponding wavelengths k* and 
m* as a function of Wi. We find that the shear-banded 
profile is unstable on the whole plateau, contrary to the 
axisymmetric prediction for our parameters [T5] . 

In order to determine whether the instability is interfa- 
cial or bulk, we employ the following method. For every 
Wi along the stress plateau, we construct an auxiliary 
Taylor- Couette system that models only the high-shear- 
rate band of the original flow. The flow in the auxiliary 
system is homogeneous (not shear-banded) and has the 
same laminar velocity and stress profiles as the original 
high-shear-rate band. This is achieved by choosing the 
gap size of the auxiliary flow cell to be equal to the width 
of the high-shear-rate band and rotating both inner and 
outer cylinders, while keeping Ri and the DJS parame- 
ters the same as above. As the result, the high-shear-rate 
band and its model version only differ in the nature of 
their outer boundary - the interface with the low-shear- 
rate band (soft) or a wall (hard), respectively. (For a 
recent discussion of the soft vs hard boundary conditions 
in the context of shear-banding, see [23].) We then com- 
pare the eigenspectra in both setups: the eigenmodes 
that appear in both systems clearly correspond to the 
bulk instability. No instability was found when the same 
procedure was applied to the low-shear-rate band. 



In Fig 3(a) we show the comparison of the shear- 



banded and (homogeneous) auxiliary systems at Wi = 29 
(near the upper end of the stress plateau) for k* and m*, 
the most unstable values of k and m. In Fig 3(c) we 



present the contour plot of the real part of the leading 
eigenvalue for different values of k and m in both se- 
tups. In spite of the different boundary conditions in the 
two systems, we observe that their leading eigenvalues 
coincide, which we attribute to the small size of the low- 
shear-rate band at high Wi and the influence of the no- 
slip boundary condition at the wall on the high-shear-rate 
band. This coincidence implies that the instability origi- 
nates in the bulk of the high-shear-rate band. Moreover, 
we find no local maxima besides the eigenvalue family 
corresponding to the most unstable bulk mode, which is 
dominant in a wide range of Wi ~ 10 — 30. The fact that 
we see a 3D instability where the axisymmetric analysis 
predicted stable flow [T5] , is perhaps not surprising since 
the first instability predicted for purely elastic Taylor- 
Couette flow of polymer solutions, which is similar to 
our bulk instability, is non-axisymmetric [24] . 

For small values of the applied shear rate the situation 
is different. As can be seen from Fig 3(b) for Wi = 3.0, 
the most unstable mode of the shear-banded profile does 
not coincide with any eigenvalue of the auxiliary system, 
and we conclude that this instability is driven by the 
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interface. We also observe, Fig |3(d)| that the shape of 
this eigenvalue family on the 5ftecr(fc,ra) contour plot is 
different from the bulk mode of Fig |3(c) , and is similar 
to what was found by Fielding for the interfacial mode 
in plane shear [T6] . 

At intermediate values of the applied shear rate, both 
the interfacial and bulk modes are unstable, as can be 



seen from Fig 3(e) The velocity profiles corresponding to 
the two peaks in the $lecr(k,m) contour plot differ funda- 
mentally: one has its maximum velocity in the vicinity of 
the interface, while the other is mostly present in the bulk 
(see Supplemental Material). The height of the bulk peak 
rapidly overtakes that of the interfacial peak as Wi is in- 
creased. The transition from the interfacial to the bulk 
mode takes place at the critical value of the Weissenberg 
number that depends on the diffusive length: Wi cr it ~ 10 
for I = 13 /im, and Wi cr it ~ 6 — 8 for I — 4 jam. 

Our use of the auxiliary Taylor- Couette systems to 
determine the nature of the instability is motivated by 
the observation that decreasing the value of the diffusion 
length I enhances both modes of instability. Indeed, this 
parameter not only controls the width of the interface, 
thus affecting the interfacial mode [16 , but also sets the 
minimal spatial extent of stress gradients in the fluid, 
and hence has a strong effect on the bulk mode as well. 
Therefore, one cannot deduce the nature of an eigenmode 
by observing how its eigenvalue changes with I. 

In summary, we presented theoretical evidence that 
the shear-banded Taylor- Couette flow is unstable with 
respect to non-axisymmetric perturbations. For param- 
eters matched to the experiments by Lerouge et al. 
[U [TUJ [11] , we find the interfacial instability only at the 
beginning of the stress plateau, while most of the plateau 
is occupied by the bulk instability. These results are con- 
sistent with the observations by Fardin et al. [11] , where 
the Taylor vortices, localised mostly in the high-shear- 
rate band, were observed in a wide range of shear rates 
on the stress plateau. Moreover, we find instabilities 
with the axial wavelengths of order of a few millimeters, 
roughly in agreement with the asymptotic (non-linear) 
wavelengths observed by Lerouge et al. [10]. Potentially, 
our prediction is further supported by recent experiments 
of Decruppe et al. [25] , who observed azimuthal undula- 
tions of the interface. However, unlike [8 , 10 , 11 , no axial 
interface perturbations were found in their experiment, 
and its relevance to our work remains an open question. 

As was noted by Fielding [15], the Taylor- Couette flow 
of worm-like micellar solutions is unique as it brings to- 
gether three types of instabilities: shear banding itself 
and the interfacial and bulk ones. The stability diagram 
proposed in [15], based on the axisymmetric analysis, in- 
cluded regions of interfacial and bulk instabilities sep- 
arated by a window of stable shear-banded flow, and 
any given micellar solution was predicted to belong to 
only one of these regions. Our results, based on the non- 
axisymmetric linear stability analysis, suggest that both 



types of instabilities can be found for a given micellar 
solution. We speculate that the whole curvature- Ni sta- 
bility diagram may be occupied by the unstable region, 
with the position of the transition between the interfacial 
and bulk modes, Wi cr it, being determined by the normal 
stresses in the high-shear-rate band: for larger values of 
Ni the transition would happen at smaller Wi cr i±. 

We would like to thank Sandra Lerouge for providing 
the rheological data used in Figjl] and Mike Cates, Marc 
Fardin and Sandra Lerouge for useful discussions. AM 
acknowledges support from the EPSRC Career Acceler- 
ation Fellowship (grant EP/I004262/1). 



SUPPLEMENTAL MATERIAL 

Position of the solution used by Lerouge et al. [10] 
on the stability diagram of Fielding [15] 

Since the relative curvature of the flow cell used by 
Lerouge et al. [10] was fixed at 1.13/13.33 « 0.08, the 
stability of the shear-banded flow as predicted by Field- 
ing [15] only depends on the model parameter 1 /(i-a). 
For this value of the curvature, the fluid was found to ex- 
hibit the interfacial instability for 1 /(i-a) < 10, the bulk 
instability for Y(i-a) > 600, or to be stable for intermedi- 
ate values of Y(i-a) [15]. The position of the boundaries 
between different types of instability are dependent on 
the parameter r] s /r] p , and the values quoted above are 
for ri s /ri P = 0.05 [T5] . 

In plane homogeneous (not banded) shear flow, the 
JS model predicts the following expression for the first 
normal stress difference Ni as a function of the applied 
shear rate 7 app : 



N, = 



2\rjp'y app 



A 2 7, 2 



app 



(i-« 2 )' 



(3) 



Eq.Q predicts that Ni approaches a horizontal asymp- 
tote as j app tends to infinity, and the normal-stress differ- 
ence of the high-shear-rate band TV-f 1 is well approximated 
by the asymptotic value of N\. In the limit a —> 1, Eq.([3| 
yields: 



X(l-a)' 



(4) 



and hence, the model parameter 1 /(i-a) is proportional to 
the value of the first normal stress difference in the high- 
shear-rate band. Using the data of S. Lerouge, Fig.l of 
the main text, we estimate N± « 3 — 5 • 10 3 Pa. The 
relaxation time in the linear regime was found to be A « 
0.23s [10 , as mentioned in the main text. From the slope 
of the ^(japp) at small applied shear rates, the total shear 
viscosity is in the range of 20 — 50Pa-s which is consistent 
with the linear rheology value 55Pa-s found in [10]. Since 
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Vs <^ Vp-> we approximate rj p by the value of the total 
zero-shear viscosity. This yields the following estimate: 

14 < — ?— < 58. (5) 
1 — a 

As mentioned above, the exact value of the boundary 
between stable and unstable regions of the stability dia- 
gram proposed by Fielding [15] depends on r] s /r]p which 
is difficult to infer experimentally. Therefore it is possible 
that the lowest end of the range ([5| belongs to the region 
of the interfacial instability, while most of the range is in 
the stable region of the stability diagram. Nevertheless, 
both are incompatible with the observation by Fardin et 
al. [IT] of a bulk instability. 

Moreover, the model parameter a = 0.985 used in this 
study [26] gives V(i-a) « 67 - well inside the stable region 
of the stability diagram of Fielding. We confirmed this 
by studying the linear stability with respect to the ax- 
isymmetric modes and found no instability. However, the 
flow does exhibit a non-axisymmetric instability, which 
is the main finding of our paper. 

Schematic view of the original and auxiliary 
Taylor- Couette setups 

In the auxiliary setup, the interface between the bands 
is replaced by a solid boundary that moves with the same 
velocity as the low-shear-rate band in the original setup, 
see Figj4] The only difference between the two setups 
is the nature of the outer boundary condition. In the 
original case, the velocities at the outer boundary of the 
high-shear-rate band have to match the corresponding 
velocities in the low-shear-rate band (a soft boundary), 
while in the auxiliary setup, the presence of the rotating 
outer wall imposes the no-slip boundary condition on the 
velocity. 




FIG. 4. The original (left) and auxiliary (right) Taylor- 
Couette setups. "L" and "H" denote the low- and high-shear- 
rate bands, respectively. 

At high Weissenberg numbers, the low-shear-rate band 
occupies only a small portion of the gap in the origi- 
nal setup, and the high-shear-rate band experiences the 
influence of the no-slip condition imposed at the outer 
wall; the interface, therefore, acts as a relatively hard 



boundary. In this case, the bulk eigenmodes in the spec- 
tra of the original and auxiliary systems coincide, as 
can be seen from Fig. 3(a). When the Weissenberg num- 
ber is decreased, corresponding to a wider low-shear-rate 
band, the boundary condition at the interface starts to 
be "softer" and the agreement between the eigenspectra 
of the original and auxiliary systems starts to deteriorate. 

Velocity profiles associated with the interfacial and 
bulk modes at Wi — 6.4 (see Fig.3(e)) 
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FIG. 5. Velocity profiles in the r — z plane for Wi — 6.4 
and / = 4 \ivn\ (top) first peak of the dispersion relation at 
k — 14.4 and m = 6, the interfacial mode; (bottom) second 
peak at k = 25 and m = 24, the bulk mode. The red lines 
indicate the position and the width of the computational do- 
main comprising the interface between the bands. 
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